Fast algorithms are essential for the solution of large systems of linear and non-linear equations that may contain millions of unknowns. Large systems may arise from the discretization of boundary integrals using the boundary element method (BEM) in physical problems governed by a Poisson or Helmholtz equation. This large category of problems includes fluid mechanics, electrostatic, and electrodynamic applications, such as problems of estimating capacitative and scattering effects.
Writing the problem in a BEM discretized manner in terms of a number N of boundary element degrees of freedom associated with the boundary elements that discretize the integral space gives a system of equations of the formAx=b,  (1)where x∈N is the vector indicating the N degrees of freedom. The N degrees of freedom are the unknown values of a first physical quantity (such as an electric or magnetic field, or a charge or current density). Every boundary element has one or more degrees of freedom assigned to it. Also, every boundary element is represented in space by one or more point locations. We refer to these point locations as discretization points or simply as points. For the rest of the text, we assume for simplicity that each boundary element has one degree of freedom and is represented in space by one discretization point, but this should not be considered limiting or restraining. The matrix A∈N×N is the BEM matrix (also called here “interaction matrix”; in other literature it is called the “kernel matrix”), and b∈N is vector of the known values of a second physical quantity (e.g. capacitance) at each of the N discretization points. The i-th and j-th entry of A define the interaction between the i-th and j-th discretization points.
The categories of solvers used for the solution of Eqn. (1) for large systems are iterative solvers and direct solvers. Iterative solvers are conjugate gradient based methods such as GMRES (the generalized minimal residual method) and Bi-CGSTAB (the biconjugate gradient stabilized method). The efficiency of iterative solvers depends on the speed of matrix-vector multiplication. The most commonly used fast matrix-vector methods include the fast multipole methods (FMM) and recursive skeletonization (RS). These methods reduce the time requirement for the matrix-vector multiplication from O(N2) to O(N). Iterative solvers have been successful and widely adopted but they have certain disadvantages: (i) The number of iterations needed to reach the solution depends on the condition number of the BEM matrix A. As a result, the solution time degrades with the condition number of A. Ill-conditioned systems are common in geometries with very close point interactions. (ii) Iterative solvers do not exploit low-rank modifications of the matrix A. Low-rank modifications may be factored using Sherman-Morison-Woodbury formula, simplifying the solution.
Direct solvers are based on BEM matrix inversion (i.e. finding A−1), instead of matrix-vector multiplication. The matrix inversion is time-wise more expensive compared to the matrix-vector multiplication, but the solution time of the linear system does not depend on the condition number of the BEM matrix. The fast matrix inversion is based on matrix compression techniques, such as recursive skeletonization. Direct solvers are advantageous compared to iterative solvers for solving systems of equations with “multiple right-hand sides”, i.e. systems of equations which, for m an integer greater than one, have the form:AX=B,  (2)where B−N×m is a matrix of m column vectors of known second physical quantities and X∈N×m is a matrix of unknown first physical quantities. Eqn. (1) is the special case of Eqn. (2) for m=1. Each column of X is referred to as a “right hand side”. In systems with multiple right hand sides the matrix inversion A−1 only has to be computed once. The inverted matrix is reapplied to every right hand side of X. Thus, the direct solver reuses more pre-computed information compared to an iterative solver that solves the m systems independently without information reuse.
The recursive skeletonization computes a compressed representation of A or A−1 by utilizing the block low-rank structure of the matrix [1]. The same block low-rank structure of the BEM matrix is utilized in the Fast Multipole Method (FMM). For a matrix to be compressible, it has to be block-separable. All off-diagonal blocks are low rank and admit to the low-rank factorization Ai,j=LiMi,jRi,j, where Ai,j∈ni×nj denotes the i-th and j-th off-diagonal block of the matrix, Li∈ni×ri, Mi,j∈ri×rj and Rj∈rj×nj for some integers ri and rj, for ri, rj<<ni. If a hierarchy is imposed on the rows and columns of A and all levels of that hierarchy are block separable, then the matrix is compressible by recursive skeletonization. The hierarchy, also known as the skeleton structure, is imposed by a spatial decomposition of the domain of the points. In recursive skeletonization, the skeleton structure is realized by an octree. The block-separability of A depends on the adaptivity of the octree to the geometry of the points. For highly irregular geometries with close point interactions the adaptivity of the octree has been proved to be inefficient, as the off-diagonal blocks of the matrix A may be almost full rank.